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Exchange splitting of the interaction energy and the multipole expansion 
of the wave function 

Piotr Gniewek^'E3 and Bogumit Jeziorski^ RI 

Faculty of Chemistry, University of Warsaw, Pasteura 1, 02-093 Warsaw, Poland 
(Dated: 9 July 2015) 

The exchange splitting J of the interaction energy of the hydrogen atom with a proton is calculated using 
the conventional surface-integral formula Tsurf[<p], the volume-integral formula of the symmetry-adapted per¬ 
turbation theory JsAPxi^p], and a variational volume-integral formula Jvar[<p]- The calculations are based 
on the multipole expansion of the wave function (p, which is divergent for any internuclear distance R. 
Nevertheless, the resulting approximations to the leading coefficient jo in the large-i? asymptotic series 
J(i?) = 2e~^~^R{jo jiR~^ j 2 R~^ -!-■••) converge, with the rate corresponding to the convergence radii 
equal to 4, 2, and 1 when the Jvar[<p], TsurfM, and Jsapt[<p] formulas are used, respectively. Additionally, we 
observe that also the higher j^, coefficients are predicted correctly when the multipole expansion is used in the 
Tvar[<p] and JsurfM formulas. The SAPT formula JsAPTi^p] predicts correctly only the first two coefficients, 
jo and ji, gives a wrong value of jj, and diverges for higher jn- Since the variational volume-integral formula 
can be easily generalized to many-electron systems and evaluated with standard basis-set techniques of quan¬ 
tum chemistry, it provides an alternative for the determination of the exchange splitting and the exchange 
contribution of the interaction potential in general. 

PACS numbers: 31.15.-p,31.10.-|-z 


I. INTRODUCTION 

Multipole expansion of the interaction energy is one of 
the cornerstones of the theory of atomic interactionsii^. 
It is indispensable for the description of potential en¬ 
ergy curves in the domain of large interatomic separa¬ 
tions R, where it approximates the interaction energy 
Eint{R) asymptotically^ii as a series in the inverse inte¬ 
ger powers of R, 

( 1 ) 

n 

the coefficients C„ being usually referred to as the van 
der Waals constants. The series on the r.h.s. of Eq. 0 
cannot converge to Eint{R), since the latter contains also 
exponentially decaying terms due to charge penetration 
and resonance tunneling (exchange) of electrons between 
interacting systems^i^. Moreover, the multipole expan¬ 
sion o is believed to be divergent, although this has 
been rigorously proved only for i.e. for the interac¬ 
tion of a hydrogen atom with a proton^^— and for the 
second order interaction energy of two hydrogen atoms^. 
One may note, however, that the multipole expansion 
has been proved to converge for interactions of confined 
atoms^ and in calculations with finite basis setsi^, al¬ 
though the limits obtained in both these cases differ from 
the true interaction energy. 

The multipole expansion o is closely related to and 
can be obtained from the multipole expansion of the wave 


function^ 

P ^ P0+^R~'^Pn, ( 2 ) 

n 

where po is the product of the wave functions of the 
non-interacting monomers, and PnS are the multipole 
corrections to the wave function^. The r.h.s. of Eq. @ 
does not represent the asymptotic approximation of the 
exact wave function if since it lacks the full permuta- 
tional and/or spacial symmetry of ^|J. It has been shown, 
however, that except for pathological cases^i, the cor¬ 
rect asymptotic expansion of the wave function can be 
obtained by symmetry projection^, i.e., 

N 

n—1 

where A is the projection operator imposing the appro¬ 
priate symmetry of the state ip. Equation shows that 
Eq. (2]) provides in fact the asymptotic expansion of a 
genuine primitive function p, as defined, e.g., in Ref. [H. 
Similarly as the expansion of Eq. o, the multipole ex¬ 
pansions for the wave function, Eq. ©, and for the prim¬ 
itive wave function, Eq. are expected to be divergent 
in the norm (the convergence of the expansion for 
the wave function would imply the convergence for the 
interaction energy). 

The interaction energies of two or more states have 
the same asymptotic expansion when they differ only by 
exponentially small exchange terms. Eor instance, for 
the ion, for the H 2 molecule and alkali dimers, or for 
homonuclear ions with one electron outside the closed 
shell, the interaction energies for the lowest gerade and 
ungerade states can be written as 

Sg^u{R) = Q{R)±J{R). 
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where Q{R) is the so-called Coulomb energy, assumed 
to be well represented at large R by the series o, and 
J{R) is the so called exchange energy falling off exponen¬ 
tially with R (the -I- and — signs are used for the g and 
u states, respectively). The exchange energy J (or the 
exchange splitting —2J) is of paramount importance for 
the understanding of weak intermolecular interactions^^, 
chemical bonds or magnetismi^, and is relevant exper¬ 
imentally, as it determines the rates of resonant charge 
exchange processes in slow atomic collision o^^ds . 

The exchange energy J{R) and its large R asymp¬ 
totic expansion are much more difficult to calculate than 
the long-range part of the interaction energy, given by 
Eq. dl]). This is due to the fact that J(i?), as a tunneling 
effect, is sensitive to the values of wave functions in the 
classically forbidden region of the configuration space, 
where the wave function amplitudes are very small and 
are hard to determine accurately using the conventional 
techniques of electronic structure theory. Instead of the 
wave functions ipg and ij^u, it is more convenient to work 
with a primitive functioni^ (/?, such that ijjg = Ag^, and 
ipu = Au^p, where Ag, and Au are appropriate symme¬ 
try projectors. When ip is known, J{R) can be obtained 
from the surface integral formul a^^d*^!^^ , 

where M indicates the median plane of the molecule, and 
the volume integral in the denominator is taken over the 
half of the space right to M {ip is understood to be local¬ 
ized on the left of M). Atomic units are used in Eq. ([S]) 
and throughout the paper. The surface-integral method 
has been extended to hydrogen moleculei^Ti^i alkali-metal 
dimer cations^^^— , neutral homo- and heterodimers22r^^ 
excited states of Hj ion^, and interactions of diatomic 
molecules with atomic ions^^. For the discussion of other 
extensions of this theory see Ref. [s^ 

In our previous worb^ we presented a volume inte¬ 
gral formula for J{R) rooted in the Symmetry Adapted 
Perturbation Theory (SAPT), 



dsAPxM 


{ipo\VPip){ipo\p) - {ipo\Vip){ipo\Pp) 

{ip,W - {ip,\Pip)^ ’ 


where V is the operator collecting Coulombic interac¬ 
tions of particles of one monomer with those of the other, 
and P is the operator inverting the electronic coordinates 
with respect to the midpoint of the internuclear axis (for 
and H 2 ) or permuting electrons between monomers 
(for larger systems). When ip is approximated via basis 
set expansions or via the expansion in powers of the in¬ 
teraction operator V, the formula of Eq. (jb)) was shown to 
give much better results^ than the surface-integral for¬ 
mula ([S]). In fact, with p expanded in powers of V, the 
formula (jS]) gives the expansion of the exchange energy in 
the symmetrized Rayleigh-Schrbdinger (SRS) perturba¬ 
tion theory^, which forms the basis for the calculation of 
exchange effects in most of the practical implementations 
of SAPT^-i^. 


In this communication we consider another volume- 
integral formula, which is variational in its origin, and, 
as we shall show, surpasses Jsurf)^?] and Jsapt[</?] in ac¬ 
curacy, 

{p\HPp){p\p) - {p\Hp){p\Pp) 

(V?|(p)2 - {p\Pp)'^ 

A similar expression and its simplified forms have been 
considered in the literatur o^^'^^i^^d^ in the theory of res¬ 
onant and non-resonant atom-ion charge exchange pro¬ 
cesses (Landau-Zener theory) but not in the present con¬ 
text of accurate ab initio calculations of the exchange 
energy. 

The purpose of the present work is to investigate the 
performance of the surface-integral formula, Eq. dSD, and 
the two volume-integral formulas, Eqs. ([6]), and 0, in the 
prediction of the exchange splitting energy J{R) when 
the primitive function p is represented by the multipole 
expansion, of Eq. ©• Since the error of the asymptotic 
series truncated after the A^th term is the order of 
it is not obvious that this series can be useful 
to correctly generate the exponentially small exchange 
terms in the energy. In fact is has been argued^ that “it 
is hopeless to try to compute the splitting by any per¬ 
turbative series expansion in 1/R”. On the other hand, 
such a procedure, if successful, would be attractive com¬ 
putationally since the expansion ([2]) is relatively easy to 
generate and Eq. m provides the simplest approxima¬ 
tion to the wave function that is asymptotically correct 
in the whole configuration space. 

The performance of different approximations of J{R) 
can be best investigated for the molecular ion. For 
this simplest, archetypal system the exchange energy 
can be computed from the exactly known asymptotic 

. « A^—47 

expansioi>»“ — : 

J(R)=2e-^-^R(jo+ji/R + j2/R^+j3/R^ + ---)- ( 8 ) 

Cizek et al. gave accurate values of the first 52 jk’s in 
Ref. d. For the interaction of two hydrogen atoms only 
the first term in the analogous expansion is know n^^i^° , 
although even the functional form of this leading term 
has been debated recently^!. For larger systems only ap¬ 
proximate form of the leading term is knowii22“— . Its 
accuracy is hard to ascertain since no reference data suf¬ 
ficiently accurate at large R are available. Therefore, we 
performed our investigation for the molecular ion and 
compared our results with the exact formula of Eq. (|8]). 

The problem considered by us was first studied by 
Tang, Toennies, and Yiu^ who evaluated the surface in¬ 
tegral ([5]) with p represented by the multipole expansion 
of the Rayleigh-Schrbdinger (RS) perturbation series in 
V (the polarization expansior^^). They were able to 
sum this series to infinite order and have shown that the 
first, jo term in Eq. ([5]) is obtained correctly in this way. 
They also obtained reasonably good approximate values 
of ji and j 2 by representing p through the second order 
in V. These authors did not consider the alternative, 
volume-integral formulas. 


dvar[^] — 
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It may be noted that for H J the exchange energy J (R) 
can be approximately obtained directly from the series 
(HD without using the wave function expansion of Eq. ©■ 
Brezin and Zinn-Justin^ have shown that the large n 
behavior of the van der Waals constants Cn in of Eq. (HD 
is related to the exchange energy via the relation 


nOC 

Cn-- R'^-^[J{R)]'^ dR. 

Jo 


(9) 


Inserting the expansion (jH]) into Eq. (ED, performing in¬ 
tegration over R, and comparing with the known 1/n 
expansion of Crir^, 


Cn = 


1 (u-H 1)! / 2 20 

2” \ n V? 



( 10 ) 


one finds the correct values of jo = —1 and ji = —1/2, 
while for j 2 an incorrect value of 19/8 is obtained, 24% 
smaller than the accurate value of j 2 equal to 25/8. Thus 
only the first two terms in the expansion ([5]) can be ob¬ 
tained exactly in this way. This method requires the 
knowledge of an analytic form of the n dependence of 
Cm which can hardly be expected to be available for 
larger systems. 

The organization of the paper is as follows. In Sec. [TTl 
we define the multipole and polarization expansions of 
the primitive function, discuss their applications to the 
evaluation of the asymptotics of the exchange energy and 
derive analytically the series representation for jo pre¬ 
dicted by the SAPT’s volume-integral formula. In Sec. lIIIl 
we present numerical results obtained with the large- 
order multipole expansion of the primitive function and 
with the first-order polarization wave function. A suc¬ 
cessful application of a simple variational approximation 
to the primitive function is also presented. Finally, in 
Sec.HV] we present conclusions of our investigation. 


II. THEORY 

A. Primitive function and the variational volume-integral 
formula 

A primitive function is a nontrivial linear com¬ 

bination of the asymptotically degenerate wave functions 
ijjg and ipu of the gerade and ungerade states, 

(f = Ciilg + C2i>u, ( 11 ) 

from which these exact states can be recovered by pro¬ 
jection 

i’g = lO + P)ip, = 1(1 - P)(/3. (12) 

We shall also require that is a “genuine primitive 
function”—, i.e., that it is localized in the same way as 
ipo- This means that {ip\Pip) vanishes exponentially or, 
more rigorously, that 

(13) 


for all integers n > 0. 

Substituting the wave functions of the form given by 
Eq. (1121) in the Rayleigh-Ritz functional and taking into 
account that [H, P] =0 and (I±P)^ ^ I±P one finds 

p p ^ {‘P\Hv’) + {¥>\HP(p) {<f\Hip)-{(p\HP<f) 


When the subtraction in Eq. HI is carried out analyti¬ 
cally, the R independent and the long-range terms 
in the numerator cancel out in view of Eq. HSD and one 
obtains Eq. HD for the exponentially vanishing energy 
splitting J(P). The explicit form of the non-relativistic 
Hamiltonian H of the hydrogen atom at point a inter¬ 
acting with a proton at point b is H = Hq + V, where 


Ho = 



1 

Ta 



1 

P’ 


(15) 


Tq and Tf) denoting the electron-nucleus distances. The 
operator V of Eq. HID appears in the SAPT formula of 
Eq. ([6]). 


B. Multipole expansion of the primitive function 

The multipole expansion of the primitive function <p, 
Eq. ([2), is obtained if the operator V is represented as 
the sum of charge-multipole interactions 

V - p-2 V2 + R-^ V3 + R-'^V4 + (16) 
Vn = -r2~^Pn-licOs9a), (17) 

where Pn{x) is the Legendre polynomial and 9a is the 
polar angle at nucleus a. The expansion of the eigen¬ 
functions of Ho -I- E in powers of R~^ leads then to the 
following recurrence equations for the van der Waals con¬ 
stants Cn and multipole corrections to the wave function 

‘fin- 

n 

Cn — ^ ) (y^O I m) , (1^) 

m—2 

n 

{Hq Eo'jcpn = ^ ) {Cni En)V^n—m, (l^) 

m—2 

where Eg = —^ and (po = are the ground- 

state energy and the wave function of the unperturbed 
hydrogen atom. To specify (pn uniquely we assumed 
here the intermediate normalization condition of </?, i.e., 
{po\p) = 1 or, equivalently, {po\pn) = 0 for n > 0. 

One can show by induction that pn is a product of 
Po and a polynomial in and cos 0a- The functions 
P 2 and Po consist of only one partial wave each, p and 
d, respectively, while the higher ones contain more than 
one partial-wave component. We calculated the van der 
Waals constants Cn and wave function corrections pn 
up to n=150 using a computer algebra program. We 
employed the exact representations of rational numbers, 
so the calculated Cn and pn are exact. 


{p\Pp) = o{R--) 
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C. Polarization expansion of the primitive function 


One can easily show that 


Another useful approximation to the primitive function 
is obtained by a finite order polarization expansion^ (or 
polarization approximation). This is an application of 
the standard Rayleigh-Schrbdinger perturbation theory 
to the Hamiltonian partitioning H — Hq + V, resulting 
in the expansion of ip in powers of V, 

OO 

= ( 20 ) 

k=0 

The individual terms in this expansion can be obtained 
from the equation 


[n/2] 

^n=Y^ (25) 

k^l 

where [r] denotes the entire value of r. Eq. (Ull) may 
be viewed as the expansion of ipn in powers of i.e., 
the polarization expansion of ipn- We shall also need the 
multipole expansion of truncated after the R~^ 

term 

K N 

(26) 

k—0 n=0 


k 

{Ho - (21) 

i=i 

where = {tpo\V, and the recursive process is 
initiated assuming = ipo and E^^'> = Eq. To fully 
specify we also assume the intermediate normaliza¬ 
tion of if, i.e., = 0 for fc > 0. For only the 

first-order term, in the series (1^01) is known in the 
closed form^. 

For systems like or H 2 the polarization expansion 
has been demonstrated numerically to converge to the 
wave function of the ground, gerade state^i^, although 
for large R the convergence radius is only marginally 
greater than unity^iiSi^ and the convergence rate be¬ 
comes prohibitively poor. Thus, at infinite order the 
sum of the series (1201) does not represent a genuine prim¬ 
itive function since it does not satisfy neither the second 
Eq. (fT^ nor the locality condition of Eq. (fUl) . Neverthe¬ 
less the polarization expansion (EOl) provides an asymp¬ 
totic approximation of the primitive function in the sense 
that^i^ 

+ 0(R-'^(^+i)), (22) 

where is the sum of the first K+1 terms in Eq. (Uni), 
^ is a suitable symmetry projector, (1 ± P)/2 in the 
case of H^, and k = 2, when at least one of interacting 
subsystems is charged and k = 3 otherwise. 

Each term in Eq. ((20l) can be represented by the mul¬ 
tipole expansion in powers of R~^, 


Note that for K > N/2 we have = $jv, where 

N 

^n = Y (27) 

n—0 

is the the multipole expansion of Eq. m truncated after 
the R~^ term. 


D. SAPT expansion of jo 

In this subsection we shall use the SAPT formula 
"/saptM and the multipole expansion for ip, Eq. ([21), to 
analytically evaluate the leading jo term in Eq. ([8l) . It is 
easy to see that in evaluating jo, the second term in the 
numerator of Eq. ([H]) can be neglected and the denomi¬ 
nator can be replaced by 1. Thus, in view of Eq. jo 
can be obtained by analyzing the large R asymptotics 
of the expression {ipo\VP^N) or its individual terms 
R~^{(po\VP(pn)- The function ipn can be written as a 
product of (fo and a polynomial /„ in Va and cos da, i-e., 
ipn = fnipo, where /o = I, /i = 0 and 

n—1 n 

fn{ra, cos 6) = E E Pi {cos 6a) (28) 

l—O m—0 

for n > 2. Employing the prolate spheroidal coordinates 
coordinates = (^a -I- rb)/R , rj = {va — ri,)/R and inte¬ 
grating by parts one finds that 


n 


(23) {g^olVPipofn) = \Re-^ 


n(l, r])dri 


L -'-i 


The series on the r.h.s. of (1^51) cannot converge in the 
whole configuration space but is expected to provide the 
large-i? asymptotic expansion of in the norm. 
Dalgarno and Lewis^ have shown that 


where 



dv{^,v) 

d( 


dr] + 0{R ^) , 


(29) 


°° / „n—l n \ 

Vi?-" ^ + ^ g,oPn-i{cos0a). (24) 

(k) 

For fc > I the coefficients (ph ', similarly as , are prod¬ 
ucts of polynomials in and cos 9a and the function ipo- 


= (C + ?7)(V’7-2)/„(^^%^,^). (30) 

Eqs. (HH) and (|5(I1) show that the leading asymptotic term 
of TsaptI/tiT’o] depends only on the values of /„ for ^ = 1, 
i.e. on the values of /„ on the line joining the nuclei a 










5 


and b. On this line df, = 0 and Pi{cos9b) = 1- Thus, 
the sought-after asymptotics in unchanged if the function 
fn{Ta, cos0) is replaced by 

n 

/n(ra, 1 ) = X! 

m—0 


valid for arbitrary sequences Pn and qn ■ Setting 

2 

Pn = dm Qn — T ' pw j (^O) 

(n -I- l)(n -I- 2) 

and summing both sides of Eq. (IMD over n from n = 0 
up to n = iV we find 


where 


n—1 


= J2d 


•nml ■ 


Z=0 


(32) 


By inspecting Eqs. (1^ and (|5(I)) one can also see that 
the leading contribution to {pq\VP p^fn) comes from the 
highest, m = n, power in Eq. m- Thus, to evaluate the 
leading asymptotics of JsAPxl/nV^o] one can replace the 
complete expression for /„ by dnx", where dnn is denoted 
by dn for brevity. 

Inserting Eq. (l28l) into Eq. m, setting 9 = 0 and com¬ 
paring coefficients at the highest power of one obtains 
the following recurrence relation 


ndji — do T T *^2 T ■ ■ ■ T dn—2- (33) 

Using Eq. (l33l) one can easily prove that 

dn dn—1 — (dn—1 dn — 2)- (34) 

n 

Taking also into that do = 1 and di = 0 one obtains 


—4d„ _ 

^(n-f l)(n-f 2)(n-b3) 

_ 2dN+i _ 1 — 2 

" (fV + 2)(fV + 3) ^Jn + 3)!’ 


(41) 


where we took into account that do = 1, go = 1 and 
d„+i — d„ = (—l)”+^/(n -I- 1)!. When N ^ oo the first 
term on the r.h.s. vanishes and the last one has the 
limit equal to —2/e -I- 1. Thus, the series on the l.h.s. 
converges to 2/e and one finds that JsAPTi'hAr] gives the 
correct asymptotics of J{R) corresponding to jo = “1- 
Having in mind that the multipole expansion of the wave 
function, Eq. leads to the divergent expansion for 
the total interaction energy we find it quite remarkable 
that, when inserted into the SAPT formula Jsapt ['/’], this 
expansion gives the convergent series for such a subtle 
effect as the asymptotic exchange splitting {—2je)Re~^. 


III. NUMERICAL RESULTS 


dn dn.— 'X — 


(-1)' 


A. Exchange splitting from the multipole expansion of the 
(35) primitive function 


and, consequently. 


dn 


E 

m—O 


ml 


(36) 


In view of Eqs. (uni) and (1501) the leading asymptotics of 
{ipo\VPr'2<po) is given by 

(n + l)(n + 2)(n -|- 3) ’ 
(37) 

Summing up contributions from ipo up to pN we find that 
the large R asymptotics of JsAPTi'hAr] is given by 




2’^+2 


(l-b?7)^(l-ry)”dry 


f-1 


N 


E 

n=0 


_-4dn_ 

(ji l)(n 2)(77- 3) 


Re-^. 


(38) 


Since the dn coefficients converge quickly to 1/e it is clear 
that the series ((551) converges, although slowly, with the 
d’Alembert ratio |an+i/a„| equal 1 —3n“^-|-0(n“^), i.e., 
with the corresponding convergence radius equal to unity. 
To evaluate its limit one can use the obvious identity 


PniPn+l dn) — Pn+lQn+l PnQn {Pn+1 Pn^Qn+l: (39) 


Since the functions fn are polynomials in and 
cosda, the exchange splittings dsAPT[4>Ar], dsurf[4>Ar] and 
dvar[4’Ar] can be obtained in a closed form for a wide 
range of N employing computer algebra software. For 
instance, for TV = 3 we obtained 

J_,[$3] = i?e-«(-|| - + |1|^ + .. -), (42) 

TsaptI^s] = 7?e ^( - fi - ^ )’ (^^) 

J_[<I>3] = i?e-«(-|^-lfiT + ^^ + ...). (44) 

Comparing with Eq. (j8]) we see that after multiplica¬ 
tion by e/2 the successive coefficients at 1/i?^ on the 
r.h.s. of Eqs. (1551) - (|55)) represent approximations to the 
jk coefficients in the expansion ([8]). These coefficients 
computed using the function and the appropriate en¬ 
ergy expressions will be denoted by 
and j™[4>Ar]. In view of Eqs. (H^ - (1551) . for TV = 3 the 
corresponding approximations to the jo coefficient are 
—49e/144, —16e/45, and —4147e/11340, and differ from 
the exact value jo = —1 by 7.5%, 3.3%, and 0.6%, re¬ 
spectively. 
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TABLE I. Values of jo calculated from the surface- and 
volume-integral formulas with the multipole expansion of 
Eq. m truncated after the R ^ipjv term. 


N 

•SAPT 

-Jo 

•surf 

-Jo 

•var 

-Jo 

0 

0.9061 

0.6796 

0.9060 9394 

2 

0.9514 

0.8601 

0.9805 2309 

3 

0.9665 

0.9250 

0.9940 6656 

4 

0.9762 

0.9625 

0.9984 1193 

5 

0.9821 

0.9811 

0.9995 7106 

6 

0.9861 

0.9905 

0.9998 8567 

7 

0.9889 

0.9953 

0.9999 6973 

8 

0.9909 

0.9976 

0.9999 9204 

9 

0.9924 

0.9988 

0.9999 9791 

10 

0.9936 

0.9994 

0.9999 9946 

CXD 

1.0 

1.0 

1.0 


TABLE 11. Values of ji calculated 
volume-integral formulas with the 
Eq. © truncated after the R~^tpN i 

from the surface- and 
multipole expansion of 
term. 

N 

^■SAPT 

•surf 

-ji 


•var 

-ji 

0 

0.0 

0.0 


0.0 

2 

0.272 

0.573 


0.465 99 

3 

0.317 

0.595 


0.492 24 

4 

0.362 

0.632 


0.503 58 

5 

0.388 

0.612 


0.502 63 

6 

0.407 

0.585 


0.501 28 

7 

0.420 

0.559 


0.500 52 

8 

0.430 

0.539 


0.500 19 

9 

0.438 

0.525 


0.500 07 

10 

0.445 

0.515 


0.500 02 

oo 

0.5 

0.5 


0.5 

In Tables HI HU and IIIII we show how the values of 

JO; jit 

and j 2 computed from 4)Ar converge when N in- 

creases 

(the rows 

for N = 1 are 

absent since ipi = 0 

so that 4)1 = 4)o). 

For the constants jo, 

ji, and j 2 the 

fastest convergence by far is observed in 

the case of the 


variational formula, whereas the SAPT formula gives the 
slowest convergence. Moreover, the sequence 


converges to a spurious value of 55/24 Ri 2.292, instead 
of the correct one, equal to 25/8 = 3.125. We also ob¬ 
served that for A: > 3 the sequences jf[4>Ar] diverge 
while the sequences generated by the surface-integral and 
variational formulas appear to be convergent for all con¬ 
stants jk that we computed. In the case of the variational 
formula this convergence is demonstrated numerically for 
/c < 6 in Table ITVl It can be seen that the variational vol¬ 
ume formula provides excellent accuracy also for further 
terms in the expansion of Eq. ®. 

The convergence of the sequence [4 *at] , shown in 

Table m to the correct limit is proved in Sec. IIIDl while 
for the sequence the proof of convergence can 

be found in Ref. |^. A rigorous mathematical proof that 
the sequence /o“[‘hAr], shown in the last column of Ta- 
blelH converges to the correct limit is more complicated^ 
and is beyond the scope of the present communication. 

TABLE III. Values of j 2 calculated from the surface- and 
volume-integral formulas with the multipole expansion of 
Eq. © truncated after the R ^(pN term. The exact value of 

is 25/8. The sequence ['I’iv] converges to the incorrect 

value equal to 55/24 Ri 2.292. 


N 

•SAPT 

32 

•surf 

J2 

•var 

32 

0 

1.359 

0.0 

1.359 14 

2 

1.925 

1.083 

2.626 05 

3 

2.197 

2.170 

3.099 24 

4 

2.218 

2.666 

3.147 05 

5 

2.254 

3.099 

3.156 90 

6 

2.269 

3.318 

3.146 38 

7 

2.278 

3.409 

3.136 55 

8 

2.284 

3.417 

3.130 49 

9 

2.288 

3.384 

3.127 39 

10 

2.291 

3.334 

3.125 97 

OO 

2.292 

3.125 

3.125 


The convergence of the sequences presented in Tables [TTl 
and mil as well as their limits have been established nu¬ 
merically based on the calculations for very high values 
of A. 


TABLE IV. Values of the coefficients jk, k < 6, calculated from the variational formula and the multiple expansion of the 
primitive wave function, Eq. ©, truncated after n = 10, 20, and 30. Data are rounded to show two digits differing from the 
exact value. The exact values are taken from Ref. |3- 


k 


jri^2o] 

iri'l’so] 

jk 

0 

-0.999 999 46 

-0.999 999 999 999 30 

-0.999 999 999 999 999 999 20 

-1. 

1 

-0.500 022 

-0.500 000 000 13 

-0.500 000 000 000 000 35 

-0.5 

2 

3.125 97 

3.125 000 022 

3.125 000 000 000 13 

3.125 

3 

2.708 

2.729 164 0 

2.729 166 666 628 

2.729 166 666 667 

4 

10.46 

10.216 39 

10.216 145 842 

10.216 145 833 

5 

38.9 

37.847 

37.864 321 3 

37.864 322 9 

6 

65. 

114.1 

113.263 89 

113.263 65 


To characterize the convergence rates of the sequences shown in Tables Ul lIIIl we computed the increment ratios 
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PivOT'O: and PN{jT") defined as 


pn[jI^^^) 




(45) 


and similarly for PNiju^'^^), and PnU™)- 

Since is a linear functional of $Ar, the 

increment ratio PNUk^^"^) is equal to the inverse of 
the d’Alembert ratio and its 

limit when N ^ oo determines the convergence ra¬ 
dius of the series X)n The surface-integral 

and variational formulas are nonlinear functionals of $jv 
but also in these cases the increment rations PN{jk^'^^) 
and PN{j™) can be interpreted as inverses of the 
d’Alembert ratios for the series with the coefficients 

and similarly for a™. We 
shall use the convergence radii of these series, given by 
the N ^ oo limit of the corresponding increment ratios, 
to characterize the convergence of the sequences at] 

and j™''[4>Ar]- 

By least-squares fitting to the data of our numerical 
calculations we obtained the following large N represen¬ 
tation of PNiJk^’^), A: = 0,1, 2, 3: 


Pa(Jo™'')= 4 - 2iV-i+0(iV-2), 
PA(jr)= 4-10 iV-i + 0(fV-2), 
PA(j2™’')=4-18iV-i + 0(iV-2), 
PA(i3™’')=4-26iV-i + 0(7V-2), 


where the relative uncertainties of the coefficients at the 
leading and subleading terms are of the order of 10“^° or 
smaller, based on analysis of j™[4>N] for N up to 60. 

In the case of the surface integral formula we obtained 


p^(jr')=2-^2-^ + 0(4-^), 

PJv(jr')=2-4iV-i+0(7V-2), (47) 

PA(jr') = 2-8iV-i+0(7V-2), 

PA(jr^) = 2-12iV-i + 0(fV-2). 

The uncertainties of the coefficients in the above formulas 
are at most of the order of 10“^°, based on analysis of 
j™“'^[(/)Af] for N up to 150. It is noteworthy that the 
analytic form of the large N representation of PnUI'^'^^) is 
different for fc = 0 and for the higher values of k. We were 
able to verify the formulas for pnUo^'^^ ) and /9Af(j™) 
(obtained initially by fitting) by a rigorous mathematical 
derivation^. 

Eqs. (I46p and (SZl) show that the series of approxi¬ 
mations to the jk coefficients obtained by using the mul¬ 
tipole expansion in the surface-integral and variational 
volume-integral formulas converge with the convergence 
radii equal to 2 and 4, respectively. This is consistent 
with the results presented in Tables mum much faster 
convergence of the variational expression. Although the 
convergence radius of the series of approximations of jk is 
independent of k, Eqs. (gni) and (HZl) show that the rate 


of convergence deteriorates somewhat with the increase 
of /c, in agreement with the data presented in Tables Ul lI VI 
When the volume-integral formula of SAPT is used 
one obtains the following inverse d’Alembert ratios for 
the series approximations to jk- 

= I + 3 iv-i + 0(iV-2), 
PN{jf^^^) = l + 2N-^ + OiN-^), 

PN{jl^^^) = l + 2N-^ + 0{N-^), 

PnU!^^^) = l + fN-^ + 0(iV-4). 

The first of these formulas is a direct consequence of 
the explicit expression for jo"^^'^[4>Af], given in Eq. (|38l) . 
while the remaining ones were obtained by least-squares 
fitting with uncertainties amounting to at most 10“^°, 
based on analysis of for N up to 150. 

Eq. (jiSl) shows that the convergence radius of the 
SAPT expansion of the coefficients jo,- ■ ■ J 3 , is equal to 
1. In such a case the d’Alembert ratio test is inconclusive, 
but we can use the Gauss criterion^ to find out whether 
the series converges or not. According to this criterion a 
series gq + ai + a 2 + - - - with a„ > 0 and the coefficients 
behaving at large n such that 

=1-1- hn~^ + 0{n~P) (49) 

^n+l 

with p > 1, converges if h > 1 and diverges otherwise. 
We therefore can conclude that the N ^ oo limits of 
the sequences Ji^^'^[4>Ar], and do 

exist, while the sequence diverges. We also 

found that for fc > 4 the sequences are di¬ 

vergent since in this case PNUk'^^'^) is smaller than 1 at 
large N (or h < 1 if the Gauss test is applied). 

In Pigure[T]we compare the convergence of the jo term 
of the asymptotic expansion of J{R) obtained when the 
multipole expansion of the primitive function is used 
in the variational, surface-integral, or SAPT formulas. 
The regular behavior of relative errors allows for the 
use of extrapolation techniques. We used the Levin u- 
transforni^ to extrapolate the results to the N ^ oo 
limit. This method accelerates the convergence of par¬ 
tial sums Zn = ao + ai + -. - + Gn with a transformed 
sequence 




(50) 


which under certain assumptions has the same limit as 
Zri^- Relative errors of jo obtained via eight-term Levin 
M-transform are presented in Fig. [T] using lines with filled 
squares. In case of j™ and eight-term Levin u- 
transform increases the accuracy by about 4 and 3.5 or¬ 
ders of magnitude, respectively, for N > 15. On the other 
hand this method based on just eight terms seems un¬ 
able to accelerate the convergence of for N > 12. 

We found, however, that when more terms are used, a 
significant decrease of error is possible. For instance, ex¬ 
trapolation from 150 values of gives a result which 

differs from the exact value of jo = —1 by only 









at accelerating the convergence of the series investigated 
in our work, provided adequately large N is used. 



FIG. 1. (Color online) Convergence of the jo term of the 
asymptotic expansion of J{R) obtained when the multipole 
expansion of the primitive function is used in the variational, 
surface-integral, or SAPT formulas. Lines with squares indi¬ 
cate values obtained from eight-term Levin-it extrapolation. 



FIG. 2. (Color online) Accuracy of the jk coefficients in the 
asymptotic expansion of J{R) calculated using different ex¬ 
change energy formulas and the multipole expansion of the 
primitive function, Eq. summed through N=30. Lines 
with squares indicate values obtained from eight-term Levin-u 
extrapolation. 


The regular convergence with respect to N occurs not 
only for jo but also for higher coefficients, allowing for 
successful extrapolation. Figure [5] presents the decimal 
logarithms of relative errors of j’^’s calculated with the 
three exchange energy formulas investigated in our work. 
This graph shows the accuracy of the raw results obtained 
for iV = 30 and the very significant gain in accuracy due 
to the Levin extrapolation. Note that in Fig.[2]we did not 
show the results of Levin’s w-transform for j™, and 
j™, as in these cases extrapolation cannot increase accu¬ 
racy. This is because the iV = 30 terms are not sufficient 
to establish the regularity of convergence. Nevertheless, 
it can be concluded that Levin’s u-transform is efficient 


B. Exchange splitting from the multipole expansion of the 
first-order polarization function 


In practice, the simplest nontrivial approximation to 
the primitive function ip is provided by the first-order po¬ 
larization function = ipo + if^^\ defined in Sec. Ill Cl 
In the present subsection we shall find out how accurate 
values of jo can be obtained from the multipole expansion 
of ip^^\ given by Eq. (l24l) . 

The value of has been already given in the 

liter ature^. Using our notation the result of Ref. can 
be restated as follows 




2 

3 


N 


E 


4 

n(n + l){n + 2)(n -|- 3) 


(51) 


The inverse d’Alembert ratio of the series on the r.h.s. of 
Eq. (ISTl) is equal to l-|-4n“^-|-0(n“^) so its convergence 
radius is equal to 1, i.e., is the same as for the series of 
Eq. (1381) . It is disappointing that the series of Eq. (|38l) , 
obtained with the function exhibit somewhat poorer 
convergence rate that the series of Eq. (EU, obtained with 
the first-order approximation to ■ Taking into 
account that the sum of the series in Eq. dSU equals 
— 1/18, we find that = —13e/36, which dif¬ 

fers from the exact value by —1.8%. This result has been 
obtained using different method (without the multipole 
expansion) in Ref. m 

The expression for jsurf[d>|^^] can be deduced from the 
results of Tang et ah— Equation (7.25) of that work im¬ 
plies that 





(52) 


One can easily find that the increment ratio for the se¬ 
quence equals 2-1-2 N~^ + 0{N~‘^), so the cor¬ 

responding convergence radius is equal 2. The N ^ oo 
limit of the r.h.s. of Eq. dSH) is —i (ln2-I-b)^, and 
the corresponding approximate value of jo, given by 
jsurfj(j)W]^ differs by —3.3% from the exact value jo = —1. 
Thus, the first-order polarization function gives a better 
approximation to the exchange splitting when used in the 
SAPT formula than in the surface-integral formula. 

The evaluation of is somewhat more compli¬ 

cated. Using the large-i? asymptotic estimation of the 
integral 


y e '"'•rlr[Pm{cos0a)Pn{cos9b)d^r = 


(fc -l- ? + 3)! 


( 53 ) 


+ 0 ( 77 ''+'+^)] 
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one can find that 


N 


^,-varr^(l)l _ _ 2 _ 

e-^° ^ 3 ^ fc(fc + l)(fc + 2)(/c + 3) 


N 


+ 2 E 


k,l=2 


{k - l)l{l - l)l[kl{k + I + 4) - 2] 
(fc + ; + 3)! ■ 


(54) 


The double sum in Eq. (15^ and its TV —>■ oo limit can be 
worked out analytically and one obtains 


- hm 
6 N—¥oo 


989 


TT 

y 


(55) 


The resulting approximate value of jo differs from the ex¬ 
act one by only 0.12%. It may be noted that the value of 
equivalent to Eq. (1551) . has been obtained ear¬ 
lier by Chipman and Hirschfelder— without the help of 
the multipole expansion using the closed-form expression 
for 

The N convergence of the sequence on the r.h.s. of 
Eq. (l5^ turns out to be rather slow. One can show 
analytically that the ratio of the A^th and the (iV -|- l)th 
increments, defined as in Eq. (H51) with $Ar replaced by 
and denoted by behaves at large N as 


a 

>o 


IS; 

Cl. 



1 1 1 

1 

a 

r" 



- 


- 

- 


- 
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1000 


FIG. 3. 
and K 


The increment rations for K = 40 (panel a) 

= 80 (panel b) as a function of N. 


occurs at N Ki Nc, where Nc is the solution of a tran¬ 
scendental equation 


p!^^Oo“) = 1 + 6iV-i + 0(iV-2). (56) 

Thus, the sequence converges at large TV as a 

series with the convergence radius equal to 1. This seems 
to be in a disagreement with Eq. (HU), which may suggest 
a faster convergence. It turns out, however, that when 
the multipole expansion of the higher polarization func¬ 
tions is used, the iV-convergence of corre¬ 

sponds also to the convergence radius 1. We have found 
by fitting the following large-TV behavior of 

p^^Hjr) = ^ + m+4)N-^ 

- {K-l)N-\\nN)-^+0{N-\\nN)-^). 

The prefactor 2K -|- 4 in the subleading term in Eq (l57l) 
shows that although the convergence radius for each K 
is equal to 1, the rate of convergence improves with in¬ 
creasing K. This rate becomes geometric only in infinite 
order in V, i.e., at K = oo, when Eq. (l46l) holds. How¬ 
ever, we have noted in Sec. im that when 

N < 2K. One can expect, then, that 

P^nHjD = PnUD = 4 - 2 iV-1 + 0(iV-2), (58) 

in the range of N values smaller or equal to 2AT. Thus, 
i® given by Eq. (1571) at N'^2K and by Eq (1571) 
at N < 2K. This is indeed the case as shown in Fig. [31 
where the increment rations p^]^\j™) are presented for 
K=40 and K=80. The switch from the fast geometric 
convergence at low N to the slow harmonic convergence 
at high N is well seen. One can show^ that this switch 


r(TV, + l)r(TV, + 2iV + 6) ^ (2iV + 2)! 

r(2W + 3) 2^iV! ■ ^ ’ 

This equation correctly describes the behavior of , for 
instance for K = 40 and iV = 80 it gives Nc = 232.3 and 
K = 80 Nc = 507.9, respectively, in a good agreement 
with the values deduced from Fig. [31 The conclusion of 
the considerations in this subsection is that the inclu¬ 
sion of high-order effects in V is necessary to obtain a 
fast converging approximation of the exchange splitting 
energy at large R. 

C. Exchange splitting from a variational approximation to 
the primitive function 

The calculations of high-order multipole corrections 
may not be practical for many-electron systems, never¬ 
theless we believe that the methods described above can 
be adapted for larger diatomics. This can be achieved via 
the use of Rayleigh-Ritz variational calculations of the 
primitive function with appropriately localized basis set, 
i.e., with trial functions satisfying the localization condi¬ 
tion of Eq. (IT3)) . In order to test the effectiveness of this 
approach we applied the basis set employed in our previ¬ 
ous work on the same system^, but restricted to orbitals 
centered exclusively at the nucleus a. This restriction 
makes the basis set orthogonal for any R. The parameter 
fl employed in Ref. [13 to control the size of the basis set 
(defined as the maximal sum of the order of the included 
Laguerre and Legendre polynomials) is closely connected 
to the extent of the multipole expansion of (p that can 
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TABLE V. Values of jo, ji and j 2 evaluated using the surface- 
and volume-integral formulas and the primitive function ap¬ 
proximated using either the multipole expansion through the 
10th order (3>io), or the variational treatment with localized, 
fl = 10 basis 


method 

jo 

ji 

12 

Jsapt[4>io] 

-0.9935 8974 

-0.444 639 

2.2909 

dsAPTi'I’io^] 

-0.9935 8952 

-0.490 467 

1.8381 

"Aurf [*^^ 10 ] 

-0.9994 0777 

-0.515 396 

3.3341 


-0.9994 0757 

-0.538 712 

3.8320 

dvari'l’io] 

-0.9999 9946 

-0.500 022 

3.1260 

Jvar[4>r(P] 

-0.9999 9965 

-0.499 994 

3.1185 

exact 

-1.0 

-0.5 

3.125 


be recovered by the variational calculation. For example 
ipN can be recovered with the D, = N, but not with the 
n = TV — 1 basis set. We performed variational calcula¬ 
tions of for 46 internuclear distances R=60, 62,... ,150, 
followed by the application of the exchange energy formu¬ 
las Jsapt[<p]; "^surfM and JvarM and least-squares fitting 
of the obtained values of J{R) to extract the constants 
jo, ji and j 2 - The results of these calculations for Li = 10 
are given in Table |V] together with the values obtained 
with the function $ 10 , which are given for comparison. 
It is seen that the agreement of the jo values calculated 
with the multipole expansion and with the variational 
approximation to ip is excellent. In the case of the higher 
constants ji and j 2 the agreement is not so good but 
reasonable, the multipole expansion giving consistently 
better results. As expected the volume-integral formula 
of Eq. (I?!) performs best not only in the case of the mul¬ 
tipole expansion of the primitive function, but also with 
the variational approximation to this function. 

IV. SUMMARY AND CONCLUSIONS 

Tang et al4^ showed that leading constant jo in the 
asymptotic expansion of the exchange splitting J{R) = 
+ jiR~^ + j 2 R~^ + ■ ■ ■) for can be ob¬ 
tained when the surface-integral formula Jsurflv’] of Eq. 
([S]) is evaluated with the multipole expanded polarization 
series for the wave function. These authors inferred that 
the polarization series converges to the primitive func¬ 
tion ip, rather than to the fully symmetric ground state 
function ijg, but his conclusion was later shown to be in¬ 
valid, as Cwiok et. al^ gave compelling evidence of the 
convergence to 'ipg. 

Our work significantly extends the results of Tang et 
alM as we have calculated not only the leading but also 
higher terms in the asymptotic expansion of J{R). We 
applied the surface-integral formula as well, but also two 
other exchange energy expressions that have not been 
previously considered in the present context: the volume- 
integral formula JsAPTiv?], Eq. ([5]), and the powerful 


volume-integral formula Jvarlv’] based on the Rayleigh- 
Ritz variational principle, Eq. 0. We have also calcu¬ 
lated the convergence radii for the expansions resulting 
when of the constants jo,- ■ ■,jo are calculated using the 
considered exchange energy expressions and the multi¬ 
pole expansion for the primitive function ip. 

We found that the variational and surface-integral for¬ 
mulas lead to convergent expansions for the jk constants. 
The best convergence and the largest convergence radius, 
equal to 4, was found for the expansions obtained us¬ 
ing the variational volume-integral formula. The con¬ 
vergence radius corresponding to the application of the 
surface-integral formula, dsurf[v^]) was found to be equal 
to 2. In the case of the SAPT formula, we found that 
the convergence for the constants jo, ji, and j 2 is slow 
with the convergence radius equal to 1. Moreover the 
expansion for j 2 converges to an incorrect value. The ex¬ 
pansions for the higher constants jk, k > 2, generated by 
the SAPT formula turned out to be divergent. 

When the multipole expansion for the first- or any 
finite-order polarization function is used to represent the 
primitive function ip, the resulting expansion for jo has 
convergence radius equal to one. Nevertheless the rate of 
convergence improves significantly with increasing order 
of the polarization expansion. This shows the importance 
of a high-order treatment in the interaction potential V 
to obtain a good approximation for the jo constant. Our 
results do not contradict the convergence of the polariza¬ 
tion expansion to the symmetric ground-state function 
'ijjg (for which all three energy formulas are singular) but 
rely on the fact that the polarization expansion converges 
asymptotically to the primitive function p. 

Since the calculation of high-order multipole correc¬ 
tions for many-electron systems is not practical, we have 
presented an alternative method of obtaining the primi¬ 
tive function: variational calculation with appropriately 
localized basis set. We have shown that this variationally 
obtained primitive function provides excellent values of 
jo and reasonably good approximations to higher-order 
jk constants. Also in this case the variational volume 
formula provides the most accurate results. 

Our study shows that the conventional SAPT formula 
exhibits some deficiencies for the calculation of the ex¬ 
change energy at large interatomic distances of R. The 
variational volume formula leads to much faster conver¬ 
gence and significantly more accurate results when ap¬ 
plied both with the multipole expansion of the wave func¬ 
tion and with a suitable variational approximation to this 
function. Moreover this formula provides a much better 
basis set convergence of the results than the surface inte¬ 
gral formula. One can therefore conclude that the varia¬ 
tional volume-integral formula provides an attractive al¬ 
ternative for the determination of the exchange splitting 
and the exchange contribution of the interaction poten¬ 
tial in general. 
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